/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package provider.impl.genomemap.me;

import commons.core.ModelEstimator;
import commons.util.ArrayUtil;
import genomemap.data.KSA00Data;
import genomemap.model.KSA00;
import org.apache.commons.lang.ArrayUtils;

import java.util.Arrays;

/**
 * @author Susanta Tewari
 */
public class KSA00Estimator extends ModelEstimator<KSA00, KSA00Data, KSA00Data, KSA00Estimator> {

    /**
     * Constructs ...
     *
     *
     * @param d
     */
    public KSA00Estimator(final KSA00Data d) {
        super(d);
    }

    /**
     * Method description
     *
     *
     * @param data
     *
     * @return
     */
    @Override
    public KSA00Estimator factory(final KSA00Data data) {
        return new KSA00Estimator(data);
    }

    /**
     * Method description
     *
     *
     * @return
     */
    @Override
    public KSA00 getEstimate() {

        final int probe_count  = getData().getModel().getProbes().size();
        final int ch_length    = getData().getModel().getChromosomeLength();
        final int clone_length = getData().getModel().getCloneLength();
        final int clone_count  = getData().getSampleSize();
        int counter            = 0;
        final byte[][] scores  = new byte[clone_count][];

        for (final KSA00Data.Hybridization val : getData()) {
            scores[counter++] = val.getScore();
        }


        // computation starts here.
        final int[] probeOrder = ArrayUtil.createSequence(getData().getModel().getProbes().size(),
                                     false);

        /**
         * joint[i][j] = # of clones that hybridize to both i and j,
         * and joint[i][i] = 0, for all i
         * DEPENDENT PARAMS: probe_count,clone_count,data
         */
        final int[][] joint = new int[probe_count][probe_count];

        for (int probeIndex1 = 0; probeIndex1 < probe_count; probeIndex1++) {

            for (int probeIndex2 = probeIndex1 + 1; probeIndex2 < probe_count; probeIndex2++) {

                for (int cloneIndex = 0; cloneIndex < clone_count; cloneIndex++) {

                    if ((scores[cloneIndex][probeIndex1] == 1)
                            && (scores[cloneIndex][probeIndex2] == 1)) {
                        joint[probeIndex1][probeIndex2]++;
                    }
                }

                joint[probeIndex2][probeIndex1] = joint[probeIndex1][probeIndex2];
            }
        }

        final int gap = ch_length - (probe_count * clone_length);
        final double[] y = Initial_Ys(probeOrder, probe_count, clone_count, clone_length, gap,
                                      joint);

        return getData().getModel().create(Arrays.asList(ArrayUtils.toObject(y)));
    }

    /**
     * Method description
     *
     *
     * @param probeOrder
     * @param probe_count
     * @param clone_count
     * @param clone_length
     * @param gap
     * @param joint
     *
     * @return
     */
    private static double[] Initial_Ys(final int[] probeOrder, final int probe_count,
                                       final int clone_count, final int clone_length,
                                       final int gap, final int[][] joint) {

        final double[] y  = new double[probe_count + 1];
        int msum          = 0;
        int mmax          = 0;
        int mmin          = clone_count;
        int w             = 0;
        final int[] link  = new int[probe_count];
        final int[] ulink = new int[probe_count];
        int un            = 0;

        /**
         * w = number of linked probes msum = total sum of linking clones for
         * all probes mmax = Maximum # of linking clones between any 2 probes
         * mmin = Minimum # of linking clones between any 2 probes
         */
        int t1;

        for (int probeIndex = 1; probeIndex < probe_count; probeIndex++) {

            t1 = joint[probeOrder[probeIndex] - 1][probeOrder[probeIndex - 1] - 1];

            if (t1 > 0) {

                msum += t1;

                if (mmax < t1) {
                    mmax = t1;
                }

                if (mmin > t1) {
                    mmin = t1;
                }

                link[w] = probeIndex;

                w++;

            } else {

                ulink[un] = probeIndex;

                un++;
            }
        }

        /*
         * A = space for all linked probes
         */
        final double eps1 = 0.1;
        final double A    = (double) gap - ((probe_count + 1 - w) * (clone_length - eps1));

        /*
         * If A<=0, no more space is left for the linked probes. Unlinked probes
         * = gap/(nprobe-w+1). Linked probes = 0
         */
        if (A <= 0.0) {

            y[0] = gap / (double) (probe_count + 1 - w);

            for (int probeIndex = 0; probeIndex < un; probeIndex++) {
                y[ulink[probeIndex]] = y[0];
            }

            for (int probeIndex = 0; probeIndex < w; probeIndex++) {
                y[link[probeIndex]] = 0.0;
            }

        } else {


            // Else , run linear programming model
            final double eps = A / (2 * w);

            y[0] = (gap - A) / (probe_count + 1 - w);

            for (int probeIndex = 0; probeIndex < un; probeIndex++) {
                y[ulink[probeIndex]] = y[0];
            }

            if (((clone_length - eps) * w) <= A) {

                for (int probeIndex = 0; probeIndex < w; probeIndex++) {
                    y[link[probeIndex]] = clone_length - eps;
                }

            } else {

                t1 = (mmax * w) - msum;

                final int t2 = msum - (mmin * w);

                if (t1 == 0) {

                    /* which implies that t2 would be zero as well */
                    for (int probeIndex = 0; probeIndex < w; probeIndex++) {
                        y[link[probeIndex]] = A / w;
                    }
                } else {

                    final double t3 = (A - (eps * w)) / t1;
                    final double t4 = ((w * (clone_length - eps)) - A) / t2;

                    if (t3 < t4) {

                        for (int probeIndex = 0; probeIndex < w; probeIndex++) {

                            un = joint[probeOrder[link[probeIndex]] - 1][probeOrder[link[probeIndex] - 1] - 1];
                            y[link[probeIndex]] = ((A * (mmax - un)) + (eps * ((w * un) - msum)))
                                                  / t1;
                        }

                    } else {

                        for (int probeIndex = 0; probeIndex < w; probeIndex++) {

                            un = joint[probeOrder[link[probeIndex] - 1]][probeOrder[link[probeIndex] - 1] - 1];
                            y[link[probeIndex]] = (((clone_length - eps) * (msum - (w * un)))
                                                   + (A * (un - mmin))) / t2;
                        }
                    }
                }
            }
        }

        /*
         * Assigning initial value to y[n+1] in the paper
         */
        double sum = 0.0;

        for (int probeIndex = 0; probeIndex < probe_count; probeIndex++) {
            sum += y[probeIndex];
        }

        y[probe_count] = gap - sum;


        // make y`s positive:
        for (int index = 0; index < y.length; index++) {

            if (y[index] < 0.0) {
                y[index] = -y[index];
            }
        }

        return y;
    }
}
